%% TCM модуляция

% Инициализация
clear all; % удаление всех переменных, функций и ссылок
clc;       % очистка командного окна

numFig = 0; % счетчик графиков

type = 'PSK';

% Кодирование
% Описание кода для 8-PSK
trel = struct('numInputSymbols', 4, 'numOutputSymbols', 8,...
              'numStates',       4, ...
              'nextStates', [0 0 1 1; 2 2 3 3; 0 0 1 1; 2 2 3 3],...
              'outputs',    [0 4 2 6; 1 5 3 7; 2 6 0 4; 3 7 1 5]);

% Декодирование
K     = log2(trel.numStates) + 1; % кодовое ограничение
tblen = 5 * K;                    % глубина просмотра при декодировании

% Модуляция
E = 1; % энергия сигнала

q = trel.numOutputSymbols;
m = log2(q);
n = m;
k = log2(trel.numInputSymbols); % число бит информации
R = k / n;

signals = modulationSignals(type, q, E); % сигнальные точки

df = freeDistanceOfTrellis(trel, signals); % свободное расстояние
                                           % решетчатого кода

                                           
% без кодирования
qUncoding = trel.numInputSymbols;
mUncoding = log2(qUncoding);
dUncoding = 2.0 * sqrt(E) * sin(pi / qUncoding);

%% TCM

SNRdBinit = 6;   % начальное значение отношения сигнал/шум
SNRdBfin  = 10;  % конечное значение отношения сигнал/шум

nErrMax = 100; % число ошибок, при котором завершается моделирование при
               % данном значении отношения сигнал/шум. Происходит переход
               % к следующему значению отношения сигнал/шум
               
listSNRdB = [SNRdBinit, 8, 9, SNRdBfin]; % , 10, 10.25, 10.5 9, 9.5, 10

% практические вероятности ошибок при различном соотношении сигнал/шум
Pe_bit = zeros( size(listSNRdB) );

%% Инициализация графика
numFig = numFig+1;
figure(numFig);

G  = 20 * log10(df / dUncoding);

SNR       = 10.^(listSNRdB / 10);
SNRdB_bit = 10 * log10( SNR / (R * m) ); % (E/N0)_bit = (E/N0) / (R * m);

listSNRdB_Th     = SNRdBinit : 0.1 : SNRdBfin + 4;
listSNR_Th       = 10.^(listSNRdB_Th / 10);
listSNRdB_bit_Th = 10 * log10( listSNR_Th / (R * m) );
% (E/N0)_bit = (E/N0) / m;
listSNRdB_bitUncoding_Th = 10 * log10(listSNR_Th / mUncoding);

% SNR    = E/noise      => noise  = E/SNR
% DNoise = noise/2      => DNoise = E / (2*SNR)
% sigma  = sqrt(DNoise) => sigma  = sqrt( E / (2*SNR) )
twoSigma = 2 * sqrt( E ./ (2 * listSNR_Th) );

thoerPe         = Q(df ./ twoSigma);
thoerPeUncoding = 2 * Q( sqrt(2 * listSNR_Th) * sin(pi / qUncoding) );

thoerPe_bit         = thoerPe / (R * m);
thoerPe_bitUncoding = thoerPeUncoding / mUncoding;

%% Моделирование
j = 1; % счетчик отношений сигнал к шуму
for SNRdB = listSNRdB
    
    % SNR = E/No
    % SNR - signal to noise ratio (отношение сигнал/шум)
    SNR = 10^(SNRdB/10); % перевод из децибел
    
    noise  = E/SNR;        % шум
    DNoise = noise/2;      % дисперсия шума, для БГШ D = NO/2
    sigma  = sqrt(DNoise); % среднее квадратичное отклонение
    
    l = 0; % счетчик испытаний
    
    % счетчики ошибок
    nBitErr         = 0;
    nErr            = 0;
    
    finalstate = []; % состояние решетчатого кодера
    % состояния декодера Витерби
    f1 = []; f2 = []; f3 = [];
    
    msgs = zeros(k*tblen, 1); % история сообщений
    
    while nBitErr < nErrMax
        % Случайные данные
        msg = randi([0 1], [k, 1]);
        
        % Модуляция TCM
        [point, finalstate] = tcmEncoder(msg, trel, signals, finalstate);
        
        % Добавление шума
        noisePoint    = [sigma*randn; sigma*randn];
        point         = point + noisePoint;
        
        % Демодуляция TCM
        [decoded, f1, f2, f3] = tcmDecoder(point, trel, signals, tblen,...
                                           f1, f2, f3);
        
        if l >= tblen && ~isequal(msgs(1:k), decoded)
            % сосчитать число битовых ошибок
            nBitErr = nBitErr + biterr(msgs(1:k), decoded);
            nErr = nErr + 1;
            
            Pe_bit(j)         = nBitErr / ( (l - tblen + 1) * k ); % добавил + 1
            
            fprintf('\nSNRdB = %d/%d, l = %d, nBitErr = %d/%d, Pe_bit(j) = %d\n\n',...
                    SNRdB, listSNRdB(end), l, nBitErr, nErrMax, Pe_bit(j));
            

            % Вывод графика зависимости вероятности ошибки от отношения сигнал/шум
            % график зависимости вероятности ошибки от отношения сигнал/шум
            semilogy(SNRdB_bit,                Pe_bit,         'k-o',...
                listSNRdB_bitUncoding_Th, thoerPe_bitUncoding, 'k--',...
                listSNRdB_bit_Th,         thoerPe_bit,         'k-');
            title(['P_{e}^{(bit)}(E/N_{o})_{(bit)}', '; G = ', num2str(G)])
            xlabel('(E/N_{o})_{(bit)}, dB')
            ylabel('P_{e}^{(bit)}')
            legend(['TCM(', num2str(q), '-', type, ')'],...
                ['Theoretical ', num2str(qUncoding), '-', type],...
                ['Theoretical TCM(', num2str(q), '-', type, ')'], 0)
            grid on
            drawnow

                
        end
        
        % обеспечить задержку декодирования
        msgs(1:k) = msg;
        msgs      = circshift(msgs, -k); % циклически сдвинуть 
                                         % элементы массива
        
        l = l+1; % перейти к следующему испытанию
    end
    
    j = j+1; % перейти к следующему соотношению отношения сигнала к шуму
end

% Вывод графика зависимости вероятности ошибки от отношения сигнал/шум
% график зависимости вероятности ошибки от отношения сигнал/шум
semilogy(SNRdB_bit,                Pe_bit,              'k-o',...
         listSNRdB_bitUncoding_Th, thoerPe_bitUncoding, 'k--',...
         listSNRdB_bit_Th,         thoerPe_bit,         'k-');
title(['P_{e}^{(bit)}(E/N_{o})_{(bit)}', '; G = ', num2str(G)])
xlabel('(E/N_{o})_{(bit)}, dB')
ylabel('P_{e}^{(bit)}')
legend(['TCM(', num2str(q), '-', type, ')'],...
       ['Theoretical ', num2str(qUncoding), '-', type],...
       ['Theoretical TCM(', num2str(q), '-', type, ')'], 0)
grid on
drawnow
